HRS HCAP Co-Calibration
Dan Mungas
2024-10-21
Methods Overview
- Multidimensional IRT used to co-calibrate HRS and HCAP cognitive measures
- Merged HRS-HCAP data set from 2016 (HRS wave 13) used for co-calibration
- Measures evaluated
- TICS3 (wave 3)
- Word List Recall (immediate and delayed),
- Count Backwards
- Serial 7’s
- Orientation to Time
- Naming
- TICS10 (wave 3 + Number Series)
- TICS3 + Number Series
- TICS13 (in TICS variants) / TICS (Comparisons with non-TICS measures)
- TICS3 + Number Series + Animal Fluency
- TICS16
- TICS13 + Trail Making A & B
- MMSE
- HCAP
- MMSE+TICS+HCAP
- TICS3 (wave 3)
- Simulation of cognitive trajectories for different measures
- True cognition simulated based on empirical modeling of cognitive trajectories from a large sample of demographically and racially diverse older adults
- UC Davis ADRC
- Kaiser cohort studies
- KHANDLE
- STAR
- Life After 90
- Simulated item response datasets created using true cognition values and item parameters from multidimensional IRT co-calibration
- True cognition simulated based on empirical modeling of cognitive trajectories from a large sample of demographically and racially diverse older adults
Data
# ********************************** Data **************************************
hcog <- read.csv("~/Psychometrics Conference/2024/Simulation WG/Data/HRS_Merge3.csv")
names(hcog) <- tolower(names(hcog))
# MMSE
hcog <- hcog %>% mutate(
ornttime_mmse = r1orient_time,
orntplace_mmse = r1orient_place,
imrec_mmse = r1imrc3,
delrec_mmse = r1dlrc3,
spellbck_mmse = backspel,
naming_mmse = r1object,
phrase_mmse = case_when(
!is.na(r1repeat) ~ case_when(
r1repeat == 1 ~ 1,
r1repeat == 5 ~ 0
)
),
commfoll_mmse = case_when(
!is.na(r1combfol) ~ case_when(
r1combfol %in% c(1,2) ~ 1,
r1combfol == 5 ~ 0
)
),
writesent_mmse = case_when(
!is.na(r1senten) ~ case_when(
r1senten %in% c(1,2) ~ 1,
r1senten == 5 ~ 0
)
),
draw_mmse = case_when(
!is.na(r1draw) ~ case_when(
r1draw == 1 ~ 1,
r1draw == 5 ~ 0
)
),
comm3step_mmse = follow3_step,
mmsetot_mmse = r1mmse_score
)
# HRS TICS
hcog <- hcog %>% mutate(
wlimrc_tics = pd174,
wldlrc_tics = pd184,
# count_tics = pd170,
animfl_tics = pd196,
numser_tics = pd216,
serial7_tics = ser7sum,
ornttime_tics = case_when(
!is.na(pd151) ~ case_when(
pd151 == 1 ~ 1,
pd151 == 5 ~ 0
)
) +
case_when(
!is.na(pd152) ~ case_when(
pd152 == 1 ~ 1,
pd152 == 5 ~ 0
)
) +
case_when(
!is.na(pd153) ~ case_when(
pd153 == 1 ~ 1,
pd153 == 5 ~ 0
)
) +
case_when(
!is.na(pd154) ~ case_when(
pd154 == 1 ~ 1,
pd154 == 5 ~ 0
)
),
naming_tics =
case_when(
!is.na(pd155) ~ case_when(
pd155 == 1 ~ 1,
pd155 == 5 ~ 0
)
) +
case_when(
!is.na(pd156) ~ case_when(
pd156 == 1 ~ 1,
pd156 == 5 ~ 0
)
) +
case_when(
!is.na(pd157) ~ case_when(
pd157 == 1 ~ 1,
pd157 == 5 ~ 0
)
) +
case_when(
!is.na(pd158) ~ case_when(
pd158 == 1 ~ 1,
pd158 == 5 ~ 0
)
),
cntbck_tics = pd170 - ornttime_tics - naming_tics
)
# HCAP
hcog <- hcog %>% mutate(
naming_hcap =
case_when(
!is.na(r1scis) ~ case_when(
r1scis == 1 ~ 1,
r1scis == 5 ~ 0
)
) +
case_when(
!is.na(r1cactus) ~ case_when(
r1cactus == 1 ~ 1,
r1cactus == 5 ~ 0
)
) +
case_when(
!is.na(r1pres) ~ case_when(
r1pres == 1 ~ 1,
r1pres == 5 ~ 0
)
) +
case_when(
!is.na(r1elbow) ~ case_when(
r1elbow == 1 ~ 1,
r1elbow == 5 ~ 0
)
) +
case_when(
!is.na(r1hammer) ~ case_when(
r1hammer == 1 ~ 1,
r1hammer == 5 ~ 0
)
),
ceradtot_hcap = r1word_total,
ceraddel_hcap = r1word_dscore,
ceradrcg_hcap = r1wlrec_totscore,
animfl_hcap = r1verbal_score,
bmim_hcap = r1bm_immscore,
bmdel_hcap = r1bm_delscore,
lmim_hcap = r1lmb_immscore,
lmdel_hcap = r1lmb_delscore,
lmrcg_hcap = r1lmb_recoscore,
conprxim_hcap = r1cp_score,
conprxdel_hcap = r1cpdel_score,
dsmt_hcap = r1dig_score,
numser_hcap = r1ns_score,
ravens_hcap = r1rv_score,
trailsa_hcap = r1tma_score,
trailsb_hcap = r1tmb_score,
spatial_hcap =
case_when(
!is.na(r1store) ~ case_when(
r1store == 1 ~ 1,
r1store == 5 ~ 0
)
) +
case_when(
!is.na(r1point) ~ case_when(
r1point == 1 ~ 1,
r1point == 5 ~ 0
)
)
)
mmse_nms <- c("ornttime_mmse","orntplace_mmse","imrec_mmse","delrec_mmse",
"spellbck_mmse","naming_mmse","phrase_mmse","commfoll_mmse","writesent_mmse",
"draw_mmse","comm3step_mmse")
tics_nms <- c("wlimrc_tics","wldlrc_tics","animfl_tics",
"numser_tics","serial7_tics","ornttime_tics","naming_tics","cntbck_tics")
hcap_nms <- c("naming_hcap","ceradtot_hcap","ceraddel_hcap","ceradrcg_hcap",
"animfl_hcap","bmim_hcap","bmdel_hcap","lmim_hcap","lmdel_hcap","lmrcg_hcap",
"conprxim_hcap","conprxdel_hcap","dsmt_hcap","numser_hcap","ravens_hcap",
"trailsa_hcap","trailsb_hcap","spatial_hcap")
### Recode MMSE
# summary(hcog[,mmse_nms])
hcog$imrec_mmse <- missCode(hcog,"imrec_mmse",0,3)
hcog$delrec_mmse <- missCode(hcog,"delrec_mmse",0,3)
# table(hcog$ornttime_mmse)
# table(hcog$orntplace_mmse)
# table(hcog$imrec_mmse)
# table(hcog$delrec_mmse)
# table(hcog$spellbck_mmse)
# table(hcog$naming_mmse)
# table(hcog$phrase_mmse)
# table(hcog$commfoll_mmse)
# table(hcog$writesent_mmse)
# table(hcog$draw_mmse)
# table(hcog$comm3step_mmse)
hcog <- hcog %>% mutate(
ornttime_mmser = ornttime_mmse,
orntplace_mmser = orntplace_mmse,
imrec_mmser = case_when(
imrec_mmse %in% c(0,1) ~ 0,
TRUE ~ imrec_mmse - 1
),
delrec_mmser = delrec_mmse,
spellbck_mmser = spellbck_mmse,
naming_mmser = naming_mmse,
phrase_mmser = phrase_mmse,
commfoll_mmser = commfoll_mmse,
writesent_mmser = writesent_mmse,
draw_mmser = draw_mmse,
comm3step_mmser = comm3step_mmse
)
# table(hcog$imrec_mmse,hcog$imrec_mmser)
### End recode MMSE
### Recode TICS
# summary(hcog[,tics_nms])
#
# table(hcog$wlimrc_tics)
# table(hcog$wldlrc_tics)
# table(hcog$cntbck_tics)
# table(hcog$count_tics)
# table(hcog$animfl_tics)
# table(hcog$numser_tics)
# table(hcog$serial7_tics)
# table(hcog$ornttime_tics)
# table(hcog$naming_tics)
hcog <- hcog %>% mutate(
wlimrc_ticsr = case_when(
wlimrc_tics %in% c(9,10) ~ 9,
TRUE ~ wlimrc_tics
),
wldlrc_ticsr = case_when(
wldlrc_tics %in% c(9,10) ~ 9,
TRUE ~ wldlrc_tics
),
cntbck_ticsr = case_when(
cntbck_tics %in% c(0,1) ~ 0,
TRUE ~ cntbck_tics - 1
),
# count_ticsr = case_when(
# count_tics %in% c(0,1,2) ~ 0,
# TRUE ~ count_tics - 2
# ),
ornttime_ticsr = case_when(
ornttime_tics %in% c(0,1) ~ 0,
TRUE ~ ornttime_tics - 1
),
naming_ticsr = case_when(
naming_tics %in% c(0,1,2) ~ 0,
TRUE ~ naming_tics - 2
),
numser_ticsr = numser_tics,
serial7_ticsr = serial7_tics
)
# table(hcog$wlimrc_tics,hcog$wlimrc_ticsr)
# table(hcog$wldlrc_tics,hcog$wldlrc_ticsr)
# table(hcog$count_tics,hcog$count_ticsr)
# table(hcog$ornttime_tics,hcog$ornttime_ticsr)
# table(hcog$naming_tics,hcog$naming_ticsr)
hcog <- recodeOrdinal(df = hcog, varlist_orig = "animfl_tics", varlist_tr = "animfl_ticsr")
# table(hcog$animfl_tics,hcog$animfl_ticsr)
### End recode TICS
### Recode HCAP
# summary(hcog[,hcap_nms])
hcog$trailsa_hcap <- missCode(hcog,"trailsa_hcap",0,300)
hcog$trailsb_hcap <- missCode(hcog,"trailsb_hcap",0,300)
hcog$trailsa_hcapi <- -hcog$trailsa_hcap
hcog$trailsb_hcapi <- -hcog$trailsb_hcap
# summary(hcog[,c("trailsa_hcapi","trailsb_hcapi")])
#
# table(hcog$naming_hcap)
# table(hcog$ceradtot_hcap)
# table(hcog$ceraddel_hcap)
# table(hcog$ceradrcg_hcap)
# table(hcog$animfl_hcap)
# table(hcog$bmim_hcap)
# table(hcog$bmdel_hcap)
# table(hcog$lmim_hcap)
# table(hcog$lmdel_hcap)
# table(hcog$lmrcg_hcap)
# table(hcog$conprxim_hcap)
# table(hcog$conprxdel_hcap)
# table(hcog$dsmt_hcap)
# table(hcog$numser_hcap)
# table(hcog$ravens_hcap)
# table(hcog$trailsa_hcap)
# table(hcog$trailsb_hcap)
# table(hcog$spatial_hcap)
hcog <- hcog %>% mutate(
naming_hcapr = case_when(
naming_hcap %in% c(0,1,2,3) ~ 1,
TRUE ~ naming_hcap - 2
),
spatial_hcapr = spatial_hcap + 1
)
varlist_orig <- c("ceradtot_hcap","ceraddel_hcap","ceradrcg_hcap",
"animfl_hcap","bmim_hcap","bmdel_hcap","lmim_hcap","lmdel_hcap","lmrcg_hcap",
"conprxim_hcap","conprxdel_hcap","dsmt_hcap","numser_hcap","ravens_hcap",
"trailsa_hcapi","trailsb_hcapi")
varlist_tr <- paste0(varlist_orig,"r")
hcog <- recodeOrdinal(hcog, varlist_orig = varlist_orig, varlist_tr = varlist_tr)
# table(hcog$naming_hcapr)
# table(hcog$ceradtot_hcapr)
# table(hcog$ceraddel_hcapr)
# table(hcog$ceradrcg_hcapr)
# table(hcog$animfl_hcapr)
# table(hcog$bmim_hcapr)
# table(hcog$bmdel_hcapr)
# table(hcog$lmim_hcapr)
# table(hcog$lmdel_hcapr)
# table(hcog$lmrcg_hcapr)
# table(hcog$conprxim_hcapr)
# table(hcog$conprxdel_hcapr)
# table(hcog$dsmt_hcapr)
# table(hcog$numser_hcapr)
# table(hcog$ravens_hcapr)
# table(hcog$trailsa_hcapir)
# table(hcog$trailsb_hcapir)
# table(hcog$spatial_hcapr)
### End recode HCAP
saveRDS(hcog, file = "~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
write.csv(hcog, file = "~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.csv", row.names=FALSE)
# --------------------------------- End Data --------------------------------Multidimensional Item Response Theory Models
# ******************************* mirt models ********************************
hcog <- readRDS("~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
### Item names for analysis
mmse_nmsr <- c("ornttime_mmser","orntplace_mmser","imrec_mmser","delrec_mmser",
"spellbck_mmser","naming_mmser","phrase_mmser","commfoll_mmser","writesent_mmser",
"draw_mmser","comm3step_mmser")
tics_nmsr <- c("wlimrc_ticsr","wldlrc_ticsr","cntbck_ticsr","animfl_ticsr",
"numser_ticsr","serial7_ticsr","ornttime_ticsr","naming_ticsr")
hcap_nmsr <- c("naming_hcapr","ceradtot_hcapr","ceraddel_hcapr","ceradrcg_hcapr",
"animfl_hcapr","bmim_hcapr","bmdel_hcapr","lmim_hcapr","lmdel_hcapr","lmrcg_hcapr",
"conprxim_hcapr","conprxdel_hcapr","dsmt_hcapr","numser_hcapr","ravens_hcapr",
"trailsa_hcapir","trailsb_hcapir","spatial_hcapr")
hcog1 <- hcog[,c(mmse_nmsr,tics_nmsr,hcap_nmsr)]
names(hcog1)
### Unidimensional model
model_hrs_hcap_mmse <- "cog = 1-37"
res_1f <- mirt(hcog1, model_hrs_hcap_mmse,
method="MHRM",
technical=list(NCYCLES=1000))
summary(res_1f)
fit_1f <- M2(res_1f,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE)
### Multidimensional model - shared methods
model_hrs_hcap_mmse_m1 <- "cog = 1-37
recmmse = 3-4
wltics = 12-13
cerhcap = 21-23
bmhcap = 25-26
lmhcap = 27-29
cnpxhcap = 30-31
trhcap = 35-36"
res_md_1 <- mirt(hcog1, model_hrs_hcap_mmse_m1,
method="MHRM",
technical=list(NCYCLES=1000))
summary(res_md_1)
coef(res_md_1)
# par <- mod2values(res_md_1)
fit_md_1 <- M2(res_md_1,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
# constrained loadings of 2-indicator factors
res_md_1a <- mirt(hcog1, model_hrs_hcap_mmse_m1,
method="MHRM",
constrain = list(c(28,38),c(128,145),c(310,327),c(396,412),c(479,495)),
technical=list(NCYCLES=1000))
summary(res_md_1a)
coef(res_md_1a)
fit_md_1a <- M2(res_md_1a,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
### Multidimensional model - shared methods and content
model_hrs_hcap_mmse_m2 <- "cog = 1-37
recmmse = 3-4
wltics = 12-13
cerhcap = 21-23
bmhcap = 25-26
lmhcap = 27-29
cnpxhcap = 30-31
trhcap = 35-36
ornt = 1,18
nmng = 6,19,20"
res_md_2 <- mirt(hcog1, model_hrs_hcap_mmse_m2,
method="MHRM",
technical=list(NCYCLES=1000))
summary(res_md_2)
coef(res_md_2)
fit_md_2 <- M2(res_md_2,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
# constrained loadings of 2-indicator factors
res_md_2a <- mirt(hcog1, model_hrs_hcap_mmse_m2,
method="MHRM",
constrain = list(c(32,44),c(150,169),c(358,377),c(454,472),c(547,565),c(9,254)),
technical=list(NCYCLES=1000))
# par <- mod2values(res_md_2)
summary(res_md_2a)
coef(res_md_2a)
fit_md_2a <- M2(res_md_2a,calcNull = TRUE, CI = 0.90, type = "C2", na.rm = TRUE,
QMC = TRUE)
save(res_1f,res_md_1,res_md_1a,res_md_2,res_md_2a,fit_1f,fit_md_1,fit_md_1a,fit_md_2,fit_md_2a,
file = "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/co_calibration_results.Rdata")
# rm(list = ls()[grepl("^fit_",ls())])
# ------------------------------end mirt models ------------------------------Model Fit Comparisons
load("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/co_calibration_results.Rdata")
fit_1f <- fit_1f %>% mutate(model = "unidimensional") %>%
relocate(model)
fit_md_1a <- fit_md_1a %>% mutate(model = "multidimensional 1") %>%
relocate(model)
fit_md_2a <- fit_md_2a %>% mutate(model = "multidimensional 2") %>%
relocate(model)
fit_summ <- bind_rows(fit_1f, fit_md_1a, fit_md_2a)
pandoc.table(fit_summ, row.names = FALSE, caption = "Model fit comparisons")| model | M2 | df | p | RMSEA | RMSEA_5 | RMSEA_95 | SRMSR |
|---|---|---|---|---|---|---|---|
| unidimensional | 4100 | 629 | 0 | 0.06927 | 0.06723 | 0.07127 | 0.1571 |
| multidimensional 1 | 1993 | 618 | 0 | 0.04399 | 0.04182 | 0.04614 | 0.1493 |
| multidimensional 2 | 1991 | 614 | 0 | 0.04416 | 0.04199 | 0.04632 | 0.1488 |
| TLI | CFI |
|---|---|
| 0.8226 | 0.8325 |
| 0.9285 | 0.9336 |
| 0.9279 | 0.9335 |
Test Information Curves
extractMIRTParm <- function(model,max_thresh = 9) {
require(mirt)
parm <- mod2values(model)
parmw <- parm %>% dplyr::select(item, name, value) %>%
pivot_wider(id_cols = item, names_from = name, values_from = value) %>%
filter(!item == "GROUP") %>% mutate(
d1 = case_when(
is.na(d1) ~ d,
TRUE ~ d1
)
)
a <- as.matrix(parmw %>% dplyr::select(a1))
d <- as.matrix(parmw %>% dplyr::select(paste0("d",1:max_thresh)))
return(list(a,d))
}
hcog <- readRDS("~/Psychometrics Conference/2024/Simulation WG/Data/hrs_hcap_co_calibration_analytic.rds")
### Item names for analysis
mmse_nmsr <- c("ornttime_mmser","orntplace_mmser","imrec_mmser","delrec_mmser",
"spellbck_mmser","naming_mmser","phrase_mmser","commfoll_mmser","writesent_mmser",
"draw_mmser","comm3step_mmser")
tics_nmsr <- c("wlimrc_ticsr","wldlrc_ticsr","cntbck_ticsr","animfl_ticsr",
"numser_ticsr","serial7_ticsr","ornttime_ticsr","naming_ticsr")
hcap_nmsr <- c("naming_hcapr","ceradtot_hcapr","ceraddel_hcapr","ceradrcg_hcapr",
"animfl_hcapr","bmim_hcapr","bmdel_hcapr","lmim_hcapr","lmdel_hcapr","lmrcg_hcapr",
"conprxim_hcapr","conprxdel_hcapr","dsmt_hcapr","numser_hcapr","ravens_hcapr",
"trailsa_hcapir","trailsb_hcapir","spatial_hcapr")
hcog1 <- hcog[,c(mmse_nmsr,tics_nmsr,hcap_nmsr)]
mirt_mmse <- "mmse = 1-11"
mirt_tics3 <- "tics = 1-6"
mirt_tics10 <- "tics = 1-7"
mirt_tics13 <- "tics = 1-8"
mirt_tics16 <- "tics = 1-10"
mirt_hcap <- "hcap = 1-18"
mirt_all <- "cog = 1-37"
true_sim <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/simulated_longitudinal_true_cognition.rds")
item_labels <- names(hcog1)
a <- extractMIRTParm((res_md_1a))[[1]]
d <- extractMIRTParm((res_md_1a))[[2]]
df <- data.frame(simdata(a = a, d = d, N = 1000, itemtype = 'graded'))
# cat("Test Information Curve - MMSE")
mod <- mirt(data = df[,1:11], model = mirt_mmse)
info_mmse <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - MMSE")
# cat("Test Information Curve - TICS3")
mod <- mirt(data = df[,c(12:14,17:19)], model = mirt_tics3)
info_tics3 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS3")
# cat("Test Information Curve - TICS10")
mod <- mirt(data = df[,c(12:14,16:19)], model = mirt_tics10)
info_tics10 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS10")
# cat("Test Information Curve - TICS13")
mod <- mirt(data = df[,12:19], model = mirt_tics13)
info_tics13 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS13")
# cat("Test Information Curve - TICS16")
mod <- mirt(data = df[,c(12:19,35:36)], model = mirt_tics16)
info_tics16 <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - TICS16")
# cat("Test Information Curve - HCAP")
mod <- mirt(data = df[,20:37], model = mirt_hcap)
info_hcap <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - HCAP")
# cat("Test Information Curve - MMSE + TICS + HCAP")
mod <- mirt(data = df, model = mirt_all)
info_all <- plot(mod, type = 'info', ylim = c(0,40),
main = "Test Information Curve - MMSE + TICS + HCAP")
Simulation of Longitudinal Cognitive Change
Simulation Functions
require(lme4)
require(mirt)
require(tidyverse)
require(ggplot2)
# show size of environment objects
# sort( sapply(ls(),function(x){object.size(get(x))}))
# rm(list = ls()[grepl("^re",ls())])
# extractMIRTParm is a function that extracts item parameters from a mirt
# estimate multidimensional (uncorrelated dimensions) model and converts
# discrimination (a) parameters to a vector of item parameters and
# threshhold (d) parameters to a matrix of d values.
# Parameters
# model - estimated mirt model object
# max_thresh - maximum number of threshold paramters across items
# Value list with 2 elements
# a - vector of a parameters, 1 for each item
# d - matrix of d parameters, rows correspons to items, columns to threshholds
extractMIRTParm <- function(model,max_thresh = 9) {
require(mirt)
parm <- mod2values(model)
parmw <- parm %>% dplyr::select(item, name, value) %>%
pivot_wider(id_cols = item, names_from = name, values_from = value) %>%
filter(!item == "GROUP") %>% mutate(
d1 = case_when(
is.na(d1) ~ d,
TRUE ~ d1
)
)
a <- as.matrix(parmw %>% dplyr::select(a1))
d <- as.matrix(parmw %>% dplyr::select(paste0("d",1:max_thresh)))
return(list(a,d))
}
# simTraj is a function that 1) estimates a mirt model using a simulated
# item response dataset, 2) calculates factor scores based on the estimated
# model, 3) estimates linear mixed effects models with true cognition and
# simulated cognition (factor scores from simulated item responses) as dependent
# variables, time in study as a fixed effect variable, and person and
# person-by-time random effects, and 4) returns person specific estimates
# (random effects) of cognitive components slopes and intercepts.
# Parameters
# data - data frame that contains true cognition value, simulated item responses,
# id variable, and time in study variable
# sim_cog_var - label for variable within dataframe (data) for which
# trajectories are being simulated
# frml - formula specification for (lmer) longitudinal mixed effects model
# iteration = iteration number passed from simulateTrajectories
# Value - list with 3 elements:
# data - analytic dataframe (data) with longitudinal model predicted cognition
# values for each assessment
# re - data frame with estimated random effects from longitudinal model.
# Includes intercept and slope random effects for true cognition (_true)
# and simulated cognition factor scores (_sim)
# fe - data frame with estimated fixed effects from longitudinal model.
# simTraj <- function(data = df, item_labels = item_labels, mirt_mod = mirt_all,
# iteration = iteration) {
simTraj <- function(data = df, sim_cog_var = "sim_cog_all", frml = frml,
iteration = iteration) {
require(tidyverse)
require(lme4)
# # estimate mirt model using simulated responses and generate factor scores
re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
int_sim_se = NA, slope_sim_se = NA)
#
# tryCatch({
# mod <- mirt(data[,item_labels], mirt_mod)
# data$sim_cog <- fscores(mod)
# data <- data %>% relocate(sim_cog, .after = true_cog)
# }, error=function(e){return(re_empty)})
# data$itertion <- iteration
# frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
# time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
# frml <- "true_cog ~ time + agebl_75 +
# time:agebl_75 + (1 | id)"
#
tryCatch({
# mixed effects longitudinal model - true cognition
# res_long_true <- lmer(true_cog ~ time + (1 + time | id), data = data)
res_long_true <- lmer(formula = frml, data = data)
# summary(res_long_true)
# str(coef(res_long_true))
data$cog_true_pred <- predict(res_long_true, newdata = data, type = "response")
re_true <- data.frame(ranef(res_long_true))
re_true <- re_true %>%
pivot_wider(id_cols = grp, names_from = term,values_from = c(condval,condsd))
names(re_true) <- c("id","int_true","slope_true","int_true_se","slope_true_se")
re_true$id <- as.integer(levels(re_true$id))[re_true$id]
fe_true <- data.frame(t(lme4::fixef(res_long_true)))
names(fe_true) <- paste0(names(fe_true),"_true")
data$sim_cog <- data[,sim_cog_var]
#res_long_sim <- lmer(sim_cog ~ time + (1 + time | id), data = data)
frml <- sub("true_cog","sim_cog",frml)
res_long_sim <- lmer(formula = frml, data = data)
# summary(res_long_sim)
data$cog_sim_pred <- predict(res_long_sim, newdata = data, type = "response")
data <- data %>% dplyr::select(-sim_cog)
re_sim <- data.frame(ranef(res_long_sim))
re_sim <- re_sim %>%
pivot_wider(id_cols = grp, names_from = term,values_from = c(condval,condsd))
names(re_sim) <- c("id","int_sim","slope_sim","int_sim_se","slope_sim_se")
re_sim$id <- as.integer(levels(re_sim$id))[re_sim$id]
fe_sim <- data.frame(t(lme4::fixef(res_long_sim)))
names(fe_sim) <- paste0(names(fe_sim),"_sim")
re <- re_true %>% left_join(re_sim, by="id")
fe <- bind_cols(fe_true,fe_sim)
fe$iteration <- iteration
}, error=function(e){re <- re_empty})
re$iteration <- iteration
return(list(data,re,fe))
}
# simulateTrajectories is a function that 1) simulates item responses for
# true cognition vectors for different components of the HRS-HCAP cognitive test
# battery (MMSE, HRS TICS, HCAP, MMSE+TICS+HCAP), and 2) generates simulated observed (factor)
# scores for each of the cognitive components (simTraj()), 3) estimates linear mixed effects models
# with true cognition and simulated cognition (factor scores from simulated item responses)
# as dependent variables, time in study as a fixed effect variable, and person and
# person-by-time random effects (simTraj()), 4) collates person specific estimates (random effects)
# of cognitive components slopes and intercepts.
#
# Parameters
# true_sim - data frame with true cognition values, rows correspond to assessments,
# columns correspond to simulation iterations.
# a_par - vector/data frame containing a (discrimination) item parameters
# from mirt co-calibration model
# d_par - vector/data frame containing d (threshhold) item parameters
# from mirt co-calibration model
# item_labels - item labels
# iter_group - number of groups of (niter) iterations
# niter - number of simulation iterations within iteration groups
# out_dir - path to directory for output of simulation results
# frml - formula specification for (lmer) longitudinal mixed effects model
#
# Value - Saves lists of results to out_dir - there is a file for each iter_group
# that contains niter elements:
# ds - list of datasets for each niter simulations. Each dataset contains
# the true cognition value (true_cog), the simulated item responses conditional
# on true_cog, simulated item response data and cognition factor scores,
# and longitudinal mixed effect model predicted cognition values for each simulated
# assessment. Measures include MMSE, TICS, HCAP, and MMSE+TICS+HCAP.
# re_all - estimated random effects for the full HRS-HCAP cognitive test battery
# (MMSE+TICS+HCAP). Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# re_mmse - estimated random effects for the MMSE component of the HRS-HCAP
# cognitive test battery. Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# re_tics - estimated random effects for the TICS component of the HRS-HCAP
# cognitive test battery. Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# re_hcap - estimated random effects for the HCAP component of the HRS-HCAP
# cognitive test battery. Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim)
# fe_all - estimated fixed effects for mixed effect model of cognition measured
# by full HRS-HCAP cognitive test battery (MMSE+TICS+HCAP).
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
# fe_mmse - estimated fixed effects for mixed effect model of cognition measured
# by MMSE.
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
# fe_tics - estimated fixed effects for mixed effect model of cognition measured
# by TICS.
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
# fe_hcap - estimated fixed effects for mixed effect model of cognition measured
# by HCAP.
# Includes intercept and slope random effects for true
# cognition (_true) and simulated cognition factor scores (_sim).
simulateTrajectories <- function(true_sim, a_par, d_par, item_labels,
niter = 20, iter_group = 5, out_dir = "Analysis/Simulation Results/",
frml = "true_cog ~ time + (1 | id)") {
require(mirt)
require(tidyverse)
for (itgrp in 1:iter_group) {
ds <- list()
re_all <- list()
re_mmse <- list()
re_tics <- list()
re_hcap <- list()
fe_all <- list()
fe_mmse <- list()
fe_tics <- list()
fe_hcap <- list()
set.seed(092724)
for (iter in 1:niter) {
cat(paste0("Group - ",itgrp, ", Iteration - ",iter,"\n"))
# simulate item responses
iteration <- ((itgrp - 1) * niter) + iter
df <- data.frame(simdata(a = a_par, d = d_par, Theta = true_sim[,paste0("vm",iteration)],
itemtype = 'graded'))
names(df) <- item_labels
df$id <- true_sim$id
df$time <- true_sim$time
df$agebl_75 <- true_sim$agebl_75
df$true_cog <- true_sim[,paste0("vm",iteration)]
df$iteration <- iteration
df <- df %>% relocate(c(id,iteration,time,true_cog))
# true random effect intercept and slope - from calculate_re_true.R
true_re <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/true_random_effects.rds")
df <- df %>% left_join((true_re %>% dplyr::select(id,int_true_qrtl,
slope_true_qrtl)),by="id")
flag_low_response <- FALSE
for (itnm in item_labels) {
if (length(table(df[,itnm])) < 2) {
flag_low_response <- TRUE
}
}
if (flag_low_response == TRUE) {
iter = iter - 1
next
}
# mirt_models for cognition measures
mirt_mmse <- "mmse = 1-11"
mirt_tics <- "tics = 1-8"
mirt_hcap <- "hcap = 1-18"
mirt_all <- "cog = 1-37"
### estimate mirt model using simulated responses and generate factor scores
# create empty re dataframe to return if error in generation
re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
int_sim_se = NA, slope_sim_se = NA)
tryCatch({
mod <- mirt(df[,item_labels], mirt_all)
df$sim_cog_all <- fscores(mod)
mod <- mirt(df[,item_labels[1:11]], mirt_mmse)
df$sim_cog_mmse <- fscores(mod)
mod <- mirt(df[,item_labels[12:19]], mirt_tics)
df$sim_cog_tics <- fscores(mod)
mod <- mirt(df[,item_labels[20:37]], mirt_hcap)
df$sim_cog_hcap <- fscores(mod)
df <- df %>%
relocate(c(sim_cog_all,sim_cog_mmse,sim_cog_tics,sim_cog_hcap),
.after = true_cog)
}, error=function(e){return(re_empty)})
# rescale factor scores to equate metric to true cognition
df <- df %>% mutate(
sim_cog_mmse_rs = (((sim_cog_mmse - mean(sim_cog_mmse)) /
sd(sim_cog_mmse)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics_rs = (((sim_cog_tics - mean(sim_cog_tics)) /
sd(sim_cog_tics)) * sd(true_cog)) + mean(true_cog),
sim_cog_hcap_rs = (((sim_cog_hcap - mean(sim_cog_hcap)) /
sd(sim_cog_hcap)) * sd(true_cog)) + mean(true_cog),
sim_cog_all_rs = (((sim_cog_all - mean(sim_cog_all)) /
sd(sim_cog_all)) * sd(true_cog)) + mean(true_cog),
)
### calculate simulated random effects
# frml <- as.formula(frml)
# frml <- as.formula("true_cog ~ time + (1 | id)")
# frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
# time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
# MMSE+TICS+HCAP
res <- simTraj(data = df, sim_cog_var = "sim_cog_all_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_all = cog_true_pred,
cog_sim_pred_all = cog_sim_pred)
re <- res[[2]]
re$label <- "MMSE+TICS+HCAP"
re_all[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "MMSE+TICS+HCAP"
fe_all[[paste0("iteration-",iteration)]] <- fe
# MMSE
res <- simTraj(data = df, sim_cog_var = "sim_cog_mmse_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_mmse = cog_true_pred,
cog_sim_pred_mmse = cog_sim_pred)
re <- res[[2]]
re$label <- "MMSE"
re_mmse[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "MMSE"
fe_mmse[[paste0("iteration-",iteration)]] <- fe
# TICS
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics = cog_true_pred,
cog_sim_pred_tics = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS"
re_tics[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS"
fe_tics[[paste0("iteration-",iteration)]] <- fe
# HCAP
res <- simTraj(data = df, sim_cog_var = "sim_cog_hcap_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_hcap = cog_true_pred,
cog_sim_pred_hcap = cog_sim_pred)
re <- res[[2]]
re$label <- "HCAP"
re_hcap[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "HCAP"
fe_hcap[[paste0("iteration-",iteration)]] <- fe
ds[[paste0("iteration-",iteration)]] <- df
} # end for iter
saveRDS(ds, file = paste0(out_dir,"ds_iteration_group_", itgrp, ".rds"))
saveRDS(re_mmse, file = paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics, file = paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
saveRDS(re_hcap, file = paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
saveRDS(re_all, file = paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
saveRDS(fe_mmse, file = paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics, file = paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
saveRDS(fe_hcap, file = paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
saveRDS(fe_all, file = paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
} # end for itgrp
for (itgrp in 1:iter_group) {
if (itgrp == 1) {
re_mmse <- readRDS(paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
re_mmse <- re_mmse %>% bind_rows()
re_tics <- readRDS(paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
re_tics <- re_tics %>% bind_rows()
re_hcap <- readRDS(paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
re_hcap <- re_hcap %>% bind_rows()
re_all <- readRDS(paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
re_all <- re_all %>% bind_rows()
fe_mmse <- readRDS(paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
fe_mmse <- fe_mmse %>% bind_rows()
fe_tics <- readRDS(paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
fe_tics <- fe_tics %>% bind_rows()
fe_hcap <- readRDS(paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
fe_hcap <- fe_hcap %>% bind_rows()
fe_all <- readRDS(paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
fe_all <- fe_all %>% bind_rows()
} else {
re_mmse_temp <- readRDS(paste0(out_dir,"re_mmse_iteration_group_", itgrp, ".rds"))
re_mmse_temp <- re_mmse_temp %>% bind_rows()
re_mmse <- bind_rows(re_mmse, re_mmse_temp)
re_tics_temp <- readRDS(paste0(out_dir,"re_tics_iteration_group_", itgrp, ".rds"))
re_tics_temp <- re_tics_temp %>% bind_rows()
re_tics <- bind_rows(re_tics, re_tics_temp)
re_hcap_temp <- readRDS(paste0(out_dir,"re_hcap_iteration_group_", itgrp, ".rds"))
re_hcap_temp <- re_hcap_temp %>% bind_rows()
re_hcap <- bind_rows(re_hcap, re_hcap_temp)
re_all_temp <- readRDS(paste0(out_dir,"re_all_iteration_group_", itgrp, ".rds"))
re_all_temp <- re_all_temp %>% bind_rows()
re_all <- bind_rows(re_all, re_all_temp)
fe_mmse_temp <- readRDS(paste0(out_dir,"fe_mmse_iteration_group_", itgrp, ".rds"))
fe_mmse_temp <- fe_mmse_temp %>% bind_rows()
fe_mmse <- bind_rows(fe_mmse, fe_mmse_temp)
fe_tics_temp <- readRDS(paste0(out_dir,"fe_tics_iteration_group_", itgrp, ".rds"))
fe_tics_temp <- fe_tics_temp %>% bind_rows()
fe_tics <- bind_rows(fe_tics, fe_tics_temp)
fe_hcap_temp <- readRDS(paste0(out_dir,"fe_hcap_iteration_group_", itgrp, ".rds"))
fe_hcap_temp <- fe_hcap_temp %>% bind_rows()
fe_hcap <- bind_rows(fe_hcap, fe_hcap_temp)
fe_all_temp <- readRDS(paste0(out_dir,"fe_all_iteration_group_", itgrp, ".rds"))
fe_all_temp <- fe_all_temp %>% bind_rows()
fe_all <- bind_rows(fe_all, fe_all_temp)
}
}
saveRDS(re_mmse, file = paste0(out_dir,"re_mmse", ".rds"))
saveRDS(re_tics, file = paste0(out_dir,"re_tics", ".rds"))
saveRDS(re_hcap, file = paste0(out_dir,"re_hcap", ".rds"))
saveRDS(re_all, file = paste0(out_dir,"re_all", ".rds"))
saveRDS(fe_mmse, file = paste0(out_dir,"fe_mmse", ".rds"))
saveRDS(fe_tics, file = paste0(out_dir,"fe_tics", ".rds"))
saveRDS(fe_hcap, file = paste0(out_dir,"fe_hcap", ".rds"))
saveRDS(fe_all, file = paste0(out_dir,"fe_all", ".rds"))
# return(list("ds" = ds, "re_all" = re_all, "re_mmse" = re_mmse,
# "re_tics" = re_tics, "re_hcap" = re_hcap))
}
simulateTrajectories2 <- function(item_labels, niter = 20, iter_group = 5,
in_dir = "Analysis/Simulation Results/",
out_dir = "Analysis/Simulation Results/",
frml = "true_cog ~ time + (1 | id)") {
require(mirt)
require(tidyverse)
for (itgrp in 1:iter_group) {
ds <- list()
re_tics3 <- list()
re_tics10 <- list()
re_tics13 <- list()
re_tics16 <- list()
fe_tics3 <- list()
fe_tics10 <- list()
fe_tics13 <- list()
fe_tics16 <- list()
set.seed(092724)
for (iter in 1:niter) {
cat(paste0("Group - ",itgrp, ", Iteration - ",iter,"\n"))
# load simulated item responses
iteration <- ((itgrp - 1) * niter) + iter
df <- readRDS(paste0(in_dir,"ds_iteration_group_",itgrp,".rds"))[[iter]]
df <- df %>% dplyr::select(c(id:true_cog,agebl_75:slope_true_qrtl,any_of(item_labels)))
# ds_iteration_group_36.rds
#
# df <- data.frame(simdata(a = a_par, d = d_par, Theta = true_sim[,paste0("vm",iteration)],
# itemtype = 'graded'))
# names(df) <- item_labels
# df$id <- true_sim$id
# df$time <- true_sim$time
# df$agebl_75 <- true_sim$agebl_75
# df$true_cog <- true_sim[,paste0("vm",iteration)]
# df$iteration <- iteration
# df <- df %>% relocate(c(id,iteration,time,true_cog))
# # true random effect intercept and slope - from calculate_re_true.R
# true_re <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Data/true_random_effects.rds")
# df <- df %>% left_join((true_re %>% dplyr::select(id,int_true_qrtl,
# slope_true_qrtl)),by="id")
#
# flag_low_response <- FALSE
# for (itnm in item_labels) {
# if (length(table(df[,itnm])) < 2) {
# flag_low_response <- TRUE
# }
# }
# if (flag_low_response == TRUE) {
# iter = iter - 1
# next
# }
# mirt_models for cognition measures
# mirt_mmse <- "mmse = 1-11"
# mirt_tics <- "tics = 1-8"
# mirt_hcap <- "hcap = 1-18"
# mirt_all <- "cog = 1-37"
mirt_tics3 <- "tics3 = 1-6"
mirt_tics10 <- "tics10 = 1-7"
mirt_tics13 <- "tics13 = 1-8"
mirt_tics16 <- "tics16 = 1-10"
### estimate mirt model using simulated responses and generate factor scores
# create empty re dataframe to return if error in generation
re_empty <- data.frame(id = NA, int_true = NA, slope_true = NA,
int_true_se = NA,slope_true_se = NA, int_sim = NA, slope_sim = NA,
int_sim_se = NA, slope_sim_se = NA)
tryCatch({
mod <- mirt(df[,item_labels[c(1:3,6:8)]], mirt_tics3)
df$sim_cog_tics3 <- fscores(mod)
mod <- mirt(df[,item_labels[c(1:3,5:8)]], mirt_tics10)
df$sim_cog_tics10 <- fscores(mod)
mod <- mirt(df[,item_labels[1:8]], mirt_tics13)
df$sim_cog_tics13 <- fscores(mod)
mod <- mirt(df[,item_labels[1:10]], mirt_tics16)
df$sim_cog_tics16 <- fscores(mod)
df <- df %>%
relocate(c(sim_cog_tics3,sim_cog_tics10,sim_cog_tics13,sim_cog_tics16),
.after = true_cog)
}, error=function(e){return(re_empty)})
# rescale factor scores to equate metric to true cognition
df <- df %>% mutate(
sim_cog_tics3_rs = (((sim_cog_tics3 - mean(sim_cog_tics3)) /
sd(sim_cog_tics3)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics10_rs = (((sim_cog_tics10 - mean(sim_cog_tics10)) /
sd(sim_cog_tics10)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics13_rs = (((sim_cog_tics13 - mean(sim_cog_tics13)) /
sd(sim_cog_tics13)) * sd(true_cog)) + mean(true_cog),
sim_cog_tics16_rs = (((sim_cog_tics16 - mean(sim_cog_tics16)) /
sd(sim_cog_tics16)) * sd(true_cog)) + mean(true_cog),
)
### calculate simulated random effects
# frml <- as.formula(frml)
# frml <- as.formula("true_cog ~ time + (1 | id)")
# frml <- "true_cog ~ time + agebl_75 + slope_true_qrtl +
# time:agebl_75 + time:slope_true_qrtl + (1 + time | id)"
# TICS3
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics3_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics3 = cog_true_pred,
cog_sim_pred_tics3 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS3"
re_tics3[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS3"
fe_tics3[[paste0("iteration-",iteration)]] <- fe
# TICS10
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics10_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics10 = cog_true_pred,
cog_sim_pred_tics10 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS10"
re_tics10[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS10"
fe_tics10[[paste0("iteration-",iteration)]] <- fe
# TICS13
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics13_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics13 = cog_true_pred,
cog_sim_pred_tics13 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS13"
re_tics13[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS13"
fe_tics13[[paste0("iteration-",iteration)]] <- fe
# TICS16
res <- simTraj(data = df, sim_cog_var = "sim_cog_tics16_rs",frml = frml,
iteration = iteration)
df <- res[[1]]
df <- df %>% rename(cog_true_pred_tics16 = cog_true_pred,
cog_sim_pred_tics16 = cog_sim_pred)
re <- res[[2]]
re$label <- "TICS16"
re_tics16[[paste0("iteration-",iteration)]] <- re
fe <- res[[3]]
fe$label <- "TICS16"
fe_tics16[[paste0("iteration-",iteration)]] <- fe
ds[[paste0("iteration-",iteration)]] <- df
} # end for iter
saveRDS(ds, file = paste0(out_dir,"ds_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics3, file = paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics10, file = paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics13, file = paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
saveRDS(re_tics16, file = paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics3, file = paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics10, file = paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics13, file = paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
saveRDS(fe_tics16, file = paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
} # end for itgrp
for (itgrp in 1:iter_group) {
if (itgrp == 1) {
re_tics3 <- readRDS(paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
re_tics3 <- re_tics3 %>% bind_rows()
re_tics10 <- readRDS(paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
re_tics10 <- re_tics10 %>% bind_rows()
re_tics13 <- readRDS(paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
re_tics13 <- re_tics13 %>% bind_rows()
re_tics16 <- readRDS(paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
re_tics16 <- re_tics16 %>% bind_rows()
fe_tics3 <- readRDS(paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
fe_tics3 <- fe_tics3 %>% bind_rows()
fe_tics10 <- readRDS(paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
fe_tics10 <- fe_tics10 %>% bind_rows()
fe_tics13 <- readRDS(paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
fe_tics13 <- fe_tics13 %>% bind_rows()
fe_tics16 <- readRDS(paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
fe_tics16 <- fe_tics16 %>% bind_rows()
} else {
re_tics3_temp <- readRDS(paste0(out_dir,"re_tics3_iteration_group_", itgrp, ".rds"))
re_tics3_temp <- re_tics3_temp %>% bind_rows()
re_tics3 <- bind_rows(re_tics3, re_tics3_temp)
re_tics10_temp <- readRDS(paste0(out_dir,"re_tics10_iteration_group_", itgrp, ".rds"))
re_tics10_temp <- re_tics10_temp %>% bind_rows()
re_tics10 <- bind_rows(re_tics10, re_tics10_temp)
re_tics13_temp <- readRDS(paste0(out_dir,"re_tics13_iteration_group_", itgrp, ".rds"))
re_tics13_temp <- re_tics13_temp %>% bind_rows()
re_tics13 <- bind_rows(re_tics13, re_tics13_temp)
re_tics16_temp <- readRDS(paste0(out_dir,"re_tics16_iteration_group_", itgrp, ".rds"))
re_tics16_temp <- re_tics16_temp %>% bind_rows()
re_tics16 <- bind_rows(re_tics16, re_tics16_temp)
fe_tics3_temp <- readRDS(paste0(out_dir,"fe_tics3_iteration_group_", itgrp, ".rds"))
fe_tics3_temp <- fe_tics3_temp %>% bind_rows()
fe_tics3 <- bind_rows(fe_tics3, fe_tics3_temp)
fe_tics10_temp <- readRDS(paste0(out_dir,"fe_tics10_iteration_group_", itgrp, ".rds"))
fe_tics10_temp <- fe_tics10_temp %>% bind_rows()
fe_tics10 <- bind_rows(fe_tics10, fe_tics10_temp)
fe_tics13_temp <- readRDS(paste0(out_dir,"fe_tics13_iteration_group_", itgrp, ".rds"))
fe_tics13_temp <- fe_tics13_temp %>% bind_rows()
fe_tics13 <- bind_rows(fe_tics13, fe_tics13_temp)
fe_tics16_temp <- readRDS(paste0(out_dir,"fe_tics16_iteration_group_", itgrp, ".rds"))
fe_tics16_temp <- fe_tics16_temp %>% bind_rows()
fe_tics16 <- bind_rows(fe_tics16, fe_tics16_temp)
}
}
saveRDS(re_tics3, file = paste0(out_dir,"re_tics3", ".rds"))
saveRDS(re_tics10, file = paste0(out_dir,"re_tics10", ".rds"))
saveRDS(re_tics13, file = paste0(out_dir,"re_tics13", ".rds"))
saveRDS(re_tics16, file = paste0(out_dir,"re_tics16", ".rds"))
saveRDS(fe_tics3, file = paste0(out_dir,"fe_tics3", ".rds"))
saveRDS(fe_tics10, file = paste0(out_dir,"fe_tics10", ".rds"))
saveRDS(fe_tics13, file = paste0(out_dir,"fe_tics13", ".rds"))
saveRDS(fe_tics16, file = paste0(out_dir,"fe_tics16", ".rds"))
# return(list("ds" = ds, "re_all" = re_all, "re_mmse" = re_mmse,
# "re_tics" = re_tics, "re_hcap" = re_hcap))
}
# dat_cols <- c("id","time","age","agebl_75","slope_true_qrtl","int_true_qrtl")
# res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
# "cog_sim_pred_tics","cog_sim_pred_hcap")
mergeSimResults <-function(file_name = "ds_iteration_group_1.rds",
out_dir = out_dir,
dat_cols = dat_cols,
res_cols = res_cols,
iter_group = 50,
niter = 20) {
res <- (readRDS(paste0(out_dir,file_name))[[1]] %>%
dplyr::select(any_of(dat_cols)))
for (i in 1:iter_group) {
for (n in 1:niter) {
iteration <- ((i-1) * niter) + n
fl_nm <- sub("_1",paste0("_",i),file_name)
ds <- readRDS(paste0(out_dir,fl_nm))[[n]] %>%
dplyr::select(all_of(res_cols))
names(ds) <- paste0(res_cols,".",iteration)
res <- res %>% bind_cols(ds)
}
}
return(res)
}Simulated Cognition Change versus True Cognition Change Plots
# sim20 <- readRDS("~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/sim20.rds")
#
# re_all <- sim20[["re_all"]]
# re_mmse <- sim20[["re_mmse"]]
# re_tics <- sim20[["re_tics"]]
# re_hcap <- sim20[["re_hcap"]]
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/10_04_24_unconditional/"
re_all <- readRDS(paste0(out_dir,"re_all.rds"))
re_mmse <- readRDS(paste0(out_dir,"re_mmse.rds"))
re_tics <- readRDS(paste0(out_dir,"re_tics.rds"))
re_hcap <- readRDS(paste0(out_dir,"re_hcap.rds"))
re_all_summ <- bind_rows(re_all)
re_mmse_summ <- bind_rows(re_mmse)
re_tics_summ <- bind_rows(re_tics)
re_hcap_summ <- bind_rows(re_hcap)
ggplot(data=re_mmse,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - MMSE")
ggplot(data=re_tics,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - TICS")
ggplot(data=re_hcap,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - HCAP")
ggplot(data=re_all,aes(slope_true,slope_sim)) +
geom_point() +
geom_smooth() +
labs(caption = "True versus simulated cognitive change - MMSE + TICS + HCAP")
Simulated Cognitive Trajectories versus True Cognition Trajectories
# Create model predicted trajectories
# # code to create merged r1 dataset
# dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
# res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
# "cog_sim_pred_tics","cog_sim_pred_hcap")
# out_dir <- "Analysis/Simulation Results/"
#
# r1 <- mergeSimResults(file_name = "ds_iteration_group_1.rds", out_dir = out_dir,
# dat_cols = dat_cols, res_cols = res_cols, iter_group = 50, niter = 20)
# saveRDS(r1,file=paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
#
# rm(r1)
# prepare longitudinal data file
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/10_15_24/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_all","cog_sim_pred_all","cog_sim_pred_mmse",
"cog_sim_pred_tics","cog_sim_pred_hcap")
r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_all_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_all_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_mn <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_se <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics_se <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_hcap_se)) %>%
pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_hcap_mn,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_hcap_mn)) %>%
pivot_longer(cols = cog_true_pred_se:cog_sim_pred_hcap_se,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_all" ~ "MMSE+TICS+HCAP",
label == "cog_sim_pred_mmse" ~ "MMSE",
label == "cog_sim_pred_tics" ~ "TICS",
label == "cog_sim_pred_hcap" ~ "HCAP",
), levels = c("True Cognition","MMSE+TICS+HCAP","HCAP","TICS","MMSE"))
)
r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
cog_sim_pred_tics_mn,cog_sim_pred_all_mn,cog_sim_pred_tics_se,cog_sim_pred_all_se) %>%
mutate(
cog_blend_mn = case_when(
time <= 5 ~ cog_sim_pred_tics_mn,
time > 5 ~ cog_sim_pred_all_mn
),
cog_blend_se = case_when(
time <= 5 ~ cog_sim_pred_tics_se,
time > 5 ~ cog_sim_pred_all_se
)
) %>% relocate(cog_blend_mn, .after = cog_sim_pred_all_mn)
r3taa <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics_se:cog_blend_se,cog_sim_pred_all_mn)) %>%
pivot_longer(cols = c(cog_sim_pred_tics_mn,cog_blend_mn),
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3tab <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics_mn:cog_blend_mn,cog_sim_pred_all_se)) %>%
pivot_longer(cols = c(cog_sim_pred_tics_se,cog_blend_se),
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3ta <- r3taa %>% left_join((r3tab %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_sim_pred_tics" ~ "TICS",
label == "cog_blend" ~ "Blended TICS,MMSE+TICS+HCAP"
), levels = c("Blended TICS,MMSE+TICS+HCAP","TICS"))
)
r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)
# cor(r1[,6:31])
### Repeat for iteration 1:100
r1 <- r1[,1:505]
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_all_mn <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_all_se <- apply(r1[,names(r1)[grepl("sim_pred_all", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_mn <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_mmse_se <- apply(r1[,names(r1)[grepl("sim_pred_mmse", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics_se <- apply(r1[,names(r1)[grepl("sim_pred_tics", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_mn <- apply(r1[,names(r1)[grepl("sim_pred_hcap", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_hcap_se <- apply(r1[,names(r1)[grepl("sim_pred_hcap", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_hcap_se)) %>%
pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_hcap_mn,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_hcap_mn)) %>%
pivot_longer(cols = cog_true_pred_se:cog_sim_pred_hcap_se,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_all" ~ "MMSE+TICS+HCAP",
label == "cog_sim_pred_mmse" ~ "MMSE",
label == "cog_sim_pred_tics" ~ "TICS",
label == "cog_sim_pred_hcap" ~ "HCAP",
), levels = c("True Cognition","MMSE+TICS+HCAP","HCAP","TICS","MMSE"))
)
r4a <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(100)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r1,r2,r3,r3a,r3b)
### End data
pd <- position_dodge(width = 0.2) # move them .2 to the left and right
ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories (1000 Simulations) - All Cognition Measures")
ggplot(data = r4a, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories (100 Simulations) - All Cognition Measures")
r4 <- r4 %>% mutate(n_sim = 1000)
r4a <- r4a %>% mutate(n_sim = 100)
r10 <- r4 %>% bind_rows(r4a) %>% mutate(
n_sim = factor(n_sim, levels = c(100,1000))
)
ggplot(data = r10, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = n_sim)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "True versus simulated cognitive trajectories by number of simulations - All Cognition Measures")
r5 <- r4 %>% filter(cognitive_measure %in% c("MMSE+TICS+HCAP","True Cognition"))
ggplot(data = r5, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - MMSE + TICS + HCAP")
r6 <- r4 %>% filter(cognitive_measure %in% c("HCAP","True Cognition"))
ggplot(data = r6, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - HCAP")
r7 <- r4 %>% filter(cognitive_measure %in% c("TICS","True Cognition"))
ggplot(data = r7, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - TICS")
r7a <- r4a %>% filter(cognitive_measure %in% c("TICS","True Cognition"))
ggplot(data = r7a, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (100 Simulations) - TICS")
r8 <- r4 %>% filter(cognitive_measure %in% c("MMSE","True Cognition"))
ggplot(data = r8, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (1000 Simulations) - MMSE")
r8a <- r4a %>% filter(cognitive_measure %in% c("MMSE","True Cognition"))
ggplot(data = r8a, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "True versus simulated cognitive change (100 Simulations) - MMSE")
r9 <- r4 %>% filter(cognitive_measure %in% c("MMSE+TICS+HCAP","MMSE"))
ggplot(data = r9, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = "Simulated cognitive change (1000 Simulations) - MMSE + TICS + HCAP vs. MMSE")
ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = str_wrap("TICS all assessments versus TICS (assessments bl-5) blended with MMSE+TICS+HCAP (assessments 6-10) (1000 Simulations)",width = 70))
rm(r4,r4a,r5,r6,r7,r7a,r8,r8a,r9,r10)Simple versus Blended Cognitive Trajectories - TICS variants
- Simple trajectories are derived from mixed effects model predicted outcome values for a specific outcome.
- Blended trajectories are derived from mixed effects model predicted outcome values for two outcomes, with the first outcome used for an initial number of assessments (0-5) and the second outcome for subsequent assessments (6-10).
Blended Results. Blending occurs after separate mixed effects models for each outcome are estimated. Plot 1 shows model predicted trajectories for each TICS variant measure. Plot 2 show model predicted trajectories for simple TICS3 and blended TICS3-TICS16.
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/tics/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_tics3","cog_sim_pred_tics3","cog_sim_pred_tics10",
"cog_sim_pred_tics13","cog_sim_pred_tics16")
r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$cog_sim_pred_tics3_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics3", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics3_se <- apply(r1[,names(r1)[grepl("sim_pred_tics3", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics10_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics10", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics10_se <- apply(r1[,names(r1)[grepl("sim_pred_tics10", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics13_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics13", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics13_se <- apply(r1[,names(r1)[grepl("sim_pred_tics13", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$cog_sim_pred_tics16_mn <- apply(r1[,names(r1)[grepl("sim_pred_tics16", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_sim_pred_tics16_se <- apply(r1[,names(r1)[grepl("sim_pred_tics16", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:cog_sim_pred_tics16_se)) %>%
pivot_longer(cols = cog_true_pred_mn:cog_sim_pred_tics16_mn,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:cog_sim_pred_tics16_mn)) %>%
pivot_longer(cols = cog_true_pred_se:cog_sim_pred_tics16_se,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_tics3" ~ "TICS (wave 3)",
label == "cog_sim_pred_tics10" ~ "TICS (wave 3 + Number Series)",
label == "cog_sim_pred_tics13" ~ "TICS (wave 13)",
label == "cog_sim_pred_tics16" ~ "TICS (wave 16)"
), levels = c("True Cognition","TICS (wave 16)","TICS (wave 13)",
"TICS (wave 3 + Number Series)","TICS (wave 3)"))
)
r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
cog_sim_pred_tics3_mn,cog_sim_pred_tics16_mn,cog_sim_pred_tics3_se,cog_sim_pred_tics16_se) %>%
mutate(
cog_blend_mn = case_when(
time <= 5 ~ cog_sim_pred_tics3_mn,
time > 5 ~ cog_sim_pred_tics16_mn
),
cog_blend_se = case_when(
time <= 5 ~ cog_sim_pred_tics3_se,
time > 5 ~ cog_sim_pred_tics16_se
)
) %>% relocate(cog_blend_mn, .after = cog_sim_pred_tics16_mn)
r3taa <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics3_se:cog_blend_se,cog_sim_pred_tics16_mn)) %>%
pivot_longer(cols = c(cog_sim_pred_tics3_mn,cog_blend_mn),
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3tab <- r2ta %>% dplyr::select(-c(cog_sim_pred_tics3_mn:cog_blend_mn,cog_sim_pred_tics16_se)) %>%
pivot_longer(cols = c(cog_sim_pred_tics3_se,cog_blend_se),
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3ta <- r3taa %>% left_join((r3tab %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_sim_pred_tics3" ~ "TICS (wave 3)",
label == "cog_blend" ~ "Blended TICS (wave 3),TICS (wave 16)"
), levels = c("Blended TICS (wave 3),TICS (wave 16)","TICS (wave 3)"))
)
r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)
ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "Plot 1. True versus simulated cognitive trajectories (1000 Simulations) - All TICS Measures")
ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = str_wrap("Plot 2. TICS (wave 3) all assessments versus TICS (wave 3 - assessments bl-5) blended with TICS (wave 16 - assessments 6-10) (1000 Simulations)",width = 70))
Blended Data - Cognitive scores for different outcomes are blended and inputted into mixed effects models.
Blended Data. Blending of two cognitive measures into one outcome occurs before mixed effects model is estimated. Plot 1a shows model predicted trajectories for simple TICS measures and selected blended measures. Plot 2a shows model predicted trajectories for simple TICS3 and blended TICS3-TICS16.
# prepare longitudinal data file
out_dir <- "~/Psychometrics Conference/2024/Simulation WG/PsyMCA-2024-Simulation/Analysis/Simulation Results/Temp/"
dat_cols <- c("id","time","agebl_75","slope_true_qrtl","int_true_qrtl")
res_cols <- c("cog_true_pred_tics3-tics13","cog_sim_pred_tics3-tics13",
"cog_sim_pred_tics3-tics3","cog_sim_pred_tics10-tics13",
"cog_sim_pred_tics10-tics10","cog_sim_pred_tics13-tics16",
"cog_sim_pred_tics13-tics13")
r1 <- readRDS(paste0(out_dir,"model_predicted_trajectory_simulations.rds"))
r1$cog_true_pred_mn <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$cog_true_pred_se <- apply(r1[,names(r1)[grepl("cog_true", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
# summary(r1$cog_true_mn)
r1$`cog_sim_pred_tics3-tics13_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics3-tics13", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics3-tics13_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics3-tics13", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics3-tics3_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics3-tics3", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics3-tics3_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics3-tics3", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics10-tics13_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics10-tics13", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics10-tics13_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics10-tics13", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics10-tics10_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics10-tics10", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics10-tics10_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics10-tics10", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics13-tics16_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics13-tics16", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics13-tics16_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics13-tics16", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics13-tics13_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics13-tics13", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics13-tics13_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics13-tics13", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics3-tics16_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics3-tics16", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics3-tics16_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics3-tics16", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r1$`cog_sim_pred_tics16-tics16_mn` <- apply(r1[,names(r1)[grepl("sim_pred_tics16-tics16", names(r1))]], 1,
function(x) mean(x, na.rm = TRUE))
r1$`cog_sim_pred_tics16-tics16_se` <- apply(r1[,names(r1)[grepl("sim_pred_tics16-tics16", names(r1))]], 1,
function(x) sd(x, na.rm = TRUE))
r2 <- r1 %>% dplyr::select(any_of(dat_cols),contains("_mn"),contains("_se")) %>%
mutate(age = agebl_75 + 75)
# summary(r2$age)
r3a <- r2 %>% dplyr::select(-c(cog_true_pred_se:`cog_sim_pred_tics16-tics16_se`)) %>%
pivot_longer(cols = cog_true_pred_mn:`cog_sim_pred_tics16-tics16_mn`,
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3b <- r2 %>% dplyr::select(-c(cog_true_pred_mn:`cog_sim_pred_tics16-tics16_mn`)) %>%
pivot_longer(cols = cog_true_pred_se:`cog_sim_pred_tics16-tics16_se`,
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3 <- r3a %>% left_join((r3b %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_true_pred" ~ "True Cognition",
label == "cog_sim_pred_tics3-tics13" ~ "TICS3-TICS13",
label == "cog_sim_pred_tics3-tics3" ~ "TICS3",
label == "cog_sim_pred_tics10-tics13" ~ "TICS10-TICS13",
label == "cog_sim_pred_tics10-tics10" ~ "TICS10",
label == "cog_sim_pred_tics13-tics16" ~ "TICS13-TICS16",
label == "cog_sim_pred_tics13-tics13" ~ "TICS13",
label == "cog_sim_pred_tics3-tics16" ~ "TICS3-TICS16",
label == "cog_sim_pred_tics16-tics16" ~ "TICS16",
), levels = c("True Cognition","TICS3-TICS13","TICS3","TICS10-TICS13",
"TICS10","TICS13-TICS16","TICS13","TICS3-TICS16","TICS16"))
)
r4 <- r3 %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
r2ta <- r2 %>% dplyr::select(id,time,agebl_75,slope_true_qrtl,
`cog_sim_pred_tics3-tics3_mn`,`cog_sim_pred_tics3-tics16_mn`,
`cog_sim_pred_tics3-tics3_se`,`cog_sim_pred_tics3-tics16_se`) # %>%
# mutate(
# cog_blend_mn = case_when(
# time <= 5 ~ cog_sim_pred_tics3_mn,
# time > 5 ~ cog_sim_pred_tics16_mn
# ),
# cog_blend_se = case_when(
# time <= 5 ~ cog_sim_pred_tics3_se,
# time > 5 ~ cog_sim_pred_tics16_se
# )
# ) %>% relocate(cog_blend_mn, .after = cog_sim_pred_tics16_mn)
r3taa <- r2ta %>% dplyr::select(-c(`cog_sim_pred_tics3-tics3_se`,`cog_sim_pred_tics3-tics16_se`)) %>%
pivot_longer(cols = c(`cog_sim_pred_tics3-tics3_mn`,`cog_sim_pred_tics3-tics16_mn`),
names_to = "label", values_to = "cog") %>% mutate(
label = sub("_mn$","",label)
)
r3tab <- r2ta %>% dplyr::select(-c(`cog_sim_pred_tics3-tics3_mn`,`cog_sim_pred_tics3-tics16_mn`)) %>%
pivot_longer(cols = c(`cog_sim_pred_tics3-tics3_se`,`cog_sim_pred_tics3-tics16_se`),
names_to = "label", values_to = "cog_se") %>% mutate(
label = sub("_se$","",label)
)
r3ta <- r3taa %>% left_join((r3tab %>%
dplyr::select(id,time,label,cog_se)),by=c("id","time","label")) %>%
mutate(
cognitive_measure = factor(case_when(
label == "cog_sim_pred_tics3-tics3" ~ "TICS (wave 3)",
label == "cog_sim_pred_tics3-tics16" ~ "Blended TICS (wave 3),TICS (wave 16)"
), levels = c("Blended TICS (wave 3),TICS (wave 16)","TICS (wave 3)"))
)
r4ta <- r3ta %>% group_by(time,slope_true_qrtl,cognitive_measure) %>% summarise(
mean_cog = mean(cog),
se_cog = mean(cog_se) / sqrt(1000)
) %>% mutate(
lower = mean_cog - (1.96 * se_cog),
upper = mean_cog + (1.96 * se_cog),
)
rm(r2,r3,r3a,r3b,r2ta,r3ta,r3taa,r3tab)
ggplot(data = r4, aes(x = time, y = mean_cog, color = slope_true_qrtl)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
facet_wrap(vars(cognitive_measure)) +
labs(caption = "Plot 1a. True versus simulated cognitive trajectories (1000 Simulations) - All TICS Measures")
ggplot(data = r4ta, aes(x = time, y = mean_cog, color = slope_true_qrtl,
linetype = cognitive_measure)) +
geom_line() +
geom_errorbar(aes(ymin=lower, ymax=upper), width=.3, position=pd) +
labs(caption = str_wrap("Plot 2a. TICS (wave 3) all assessments versus TICS (wave 3 - assessments bl-5) blended with TICS (wave 16 - assessments 6-10) (1000 Simulations)",width = 70))